Approximate solutions and scaling transformations for quadratic solitons 
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We study quadratic solitons supported by two- and three-wave parametric interactions in x' 2 ' 
nonlinear media. Both planar and two-dimensional cases are considered. We obtain very accurate, 
'almost exact', explicit analytical solutions, matching the actual bright soliton profiles, with the help 
of a specially-developed approach, based on analysis of the scaling properties. Additionally, we use 
these approximations to describe the linear tails of solitary waves which are related to the properties 
of the soliton bound states. 



I. INTRODUCTION 



One of the rapidly expanding areas of research is the 
physics of solitons - wave packets, or self-trapped beams, 
that propagate with their profiles remaining undistorted. 
In particular, parametric solitons, composed of mutually 
trapped fundamental and harmonic waves, attract inter- 
est of researchers due to a wide range of possible applica- 
tions. In optics, for example, such solitons were observed 
in media with quadratic (or x^) nonlinearity, and their 
unique features can be utilized for the creation of all- 
optical information processing devices [jjj. In general, 
parametric solitons may form in different media which 
possess resonant quadratic nonlinearities, such as plasma, 
organic superlattices 0, Bose- Einstein condensate jj), 
etc. 

Many papers have been devoted to a theoretical anal- 
ysis of quadratic solitons j|. It was shown that bright 
solitons can be stable, and hence are of most interest 
for practical applications, on the other hand parametric 
dark solitons often exhibit modulation instability. How- 
ever, due to nontrivial features of the resonant coupling, 
there still remain some properties of the bright quadratic 
solitons which have not been thoroughly described or 
completely understood. The problem here is that the 
governing equations are not integrable, and general an- 
alytical solutions can't be constructed. Thus, the vari- 
ational method Q was widely used to find approximate 
solutions. In this approach, the free parameters control- 
ling the trial functions are found by minimizing the La- 
grangian functional. A rarely mentioned limitation here 
is that this technique imposes very strong restrictions on 
the class of the trial functions, for the parameters to be 
found in an explicit form. This becomes a real draw- 
back if these test functions are not quite suitable for the 
problem at hand (see, e.g., Ref ||). 

In this paper we introduce a different approach for ob- 
taining approximate expressions for the soliton envelopes. 
At first we choose trial functions which can precisely de- 
scribe the actual soliton profiles. This step involves the 
analysis of the scaling properties of a soliton family, i.e. 
how the envelopes are transformed as a free parameter 



(propagation constant) is altered. Secondly, a specially 
developed technique is used to find the fitting parame- 
ters. It will be demonstrated that the resulting solutions 
turn out to be extremely accurate. 

The rest of the paper is organized as following. First, 
two-wave solitons in planar structures are studied. Then, 
the analysis is extended to the case of three wave mutual 
trapping in an anisotropic medium. Finally, the proper- 
ties of two-component solitary beams propagating in a 
bulk medium are investigated. 



II. ONE DIMENSIONAL SOLITONS 



A. Two-wave solitons 



1. Basic equations and their properties 



Parametric interaction between the fundamental fre- 
quency (FF) wave and its second harmonic (SH) in the 
(1+1) dimensional case can be described by a set of cou- 
pled equations for slowly varying complex amplitudes of 
the wave packets (see also |§|§]). We consider the 
case when there is no walk-off, and then in normalized 
variables have [El : 



.du d 2 u 

% ~^~ + TT + u w = °> 

oz ox z 



dw 



dz dx 2 



(3w H — u 



(1) 



where u(x, z) and w(x, z) are the FF and SH amplitudes 
correspondingly, z is the propagation distance, and the 
parameter (3 characterizes the mismatch of the linear 
phase velocities. These equations can describe: (i) spa- 
tial beams in a slab waveguide, exhibiting diffraction in 
the transverse direction x, where a ss 2; or (ii) temporal 
pulses, where x stands for the retarded time and a > is 
the ratio of the absolute values of second-order dispersion 
coefficients (the signs are those which allow stable bright 
solutions). 
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Under certain conditions, mutual trapping of two 
waves can be achieved, when diffraction (or second-order 
dispersion) is exactly compensated for by nonlinear re- 
fraction. Such stationary propagation is observed for 
a special class of beams (or temporal wave packets) - 
solitons. To find their profiles, we look for solutions of 
Eqs. (@) in the for m 



u(x, z) — Auo 
w(x,z) = Xw ^xy/X^j e 2lXz , 



(2) 



where uq and wq are the real envelope amplitudes, and 
A > is the propagation constant. After substituting 
expressions @ into Eqs. the following system of cou- 
pled ordinary differential equations can be derived: 



dx 2 
d 2 w 
dx 2 



- u a + u w = 0, 
1 



(3) 



0. 



This solution can be improved by taking into account 
higher-order terms in a series decomposition over a small 
parameter, a -1 p5|^6| . However, this approach leads 
to cumbersome expressions which are somewhat hard to 
analyze. 

Exact analytical solutions of Eqs. (||) cannot be found 
for arbitrary values of a. Thus, in order to obtain approx- 
imate solutions for the soliton profiles, the variational ap- 
proach was used. Calculations with an ansatz in the form 
of Gaussian functions predict the power distribution be- 
tween the FF and SH components quite accurately, and 
provide a close estimation for the maximum amplitudes 
in the whole parameter range < a < oo |]J|,|l7]]. How- 
ever, the trial functions do not correspond to the actual 
wave profiles, and thus the "tails", or amplitude asymp- 
totics at x — > ±oo, are not described well. In other stud- 
ies jUI, the profiles of the trial functions are chosen as 
scaled exact solutions (||) or (fjj) with arbitrary ampli- 
tudes, but fixed relative widths for the FF and SH wave 
packets. Due to this limitation, precise results are ob- 
tained only for a ~ 1 and a — > oo respectively. 



where the only free parameter is the normalized mis- 
match a = 2a + /3/A. 

For localized waves, described by Eqs. (pi), the total 
power P and Hamiltonian H are conserved . The val- 
ues of these integral characteristics, corresponding to soli- 
tons defined by Eqs. (||), can be found as pX| ]: 

P = A 3 / 2 P , H = X 5 / 2 (H -P ). (4) 

Here the renormalized power and Hamiltonian arc Pq = 
P Uo + 2<tP Wo and H = 0.4(P Uo + aP W0 ), where 

/+oo r+oo 
u 2 dx, P Wo = / w^dx, (5) 
-oo J —oo 



2. Approximate analytical solution 

In order to construct a solution without the above men- 
tioned drawbacks, we want to take into account some 
characteristic properties of the bright solitons. Specifi- 
cally, we notice that the form of the SH envelope is the 
same at a — 1 and a — > oo [see Eqs. (|) and (@)]. More- 
over, we perform numerical simulations and observe a 
very remarkable fact: for a > 1, wq(x) remains almost 
self similar. Thus, we search for an approximate solution 
with the SH component in the form: 

w Q (x) = w m sech 2 (x/p), (8) 



We are looking for bright solitons, where the field de- 
cays at infinity. Such solutions of Eqs. (||) were found 
numerically for any a > p. ff| , |l2| ] and corresponding 
solitons were shown to be stable for mismatches a > a cr 
p3| . Here, the critical parameter value a cl - is a function 
of cr, and for spatial solitons ~ 0-2- 

The properties of this soliton family were extensively 
studied in the literature. An exact solution was found at 
a = l M: 



w (x) = uq(x)/V2 = (3/2) scch 2 (x/2) . 



(6) 



On the other hand, for large mismatches, a — > oo, the 
SH component approaches wq(x) ~ Vq(x)/ (2a). Then, 
in this so-called cascading limit, the FF wave is deter- 
mined as a solution of the nonlinear Schrodinger equation 
(NLSE) , and the wave envelopes are HQ : 



vq(x) ~ 2\fu sech(x), wq(x) ~ 2 sech 2 (x). 



(7) 



where the maximum amplitude w m and characteristic 
width p are unknown parameters. Then, the FF com- 
ponent profile can be determined using the first equation 
in (|^). This is a linear eigenvalue problem, which has an 
exact solution for the effective waveguide created by the 
SH field (||). In a single bright soliton, uo(x) doesn't have 
zeros, and thus we take the fundamental mode. This re- 
quirement leads to a relation between the parameters of 
an effective waveguide: 



1 + 1/P- 



(9) 



Then the corresponding FF profile is found to be 

u (x) = u m sech p (x/p), (10) 

where u m is the peak amplitude. 

Our trial functions do not satisfy the SH equation in 
(0) exactly, and thus it should be matched approximately. 
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This can be done with the help of the variational method. 
However, this leads to a set of transcendental equations, 
and then the solution parameters can't be expressed in 
an explicit form. As our aim is to derive a simple analyt- 
ical approximation, which should be easy to analyze and 
use in calculations, we choose another approach. First, 
in order to match the soliton peak, we require the equa- 
tion for the SH amplitude in (||) to be exactly satisfied 
at the soliton center, x = 0, and thus obtain: 



2w m (a + 2/p 2 ) = 



(11) 



Secondly, we notice that Eqs. (||) describe an equiva- 
lent dynamical problem, viz. particle movement with 
generalized velocities (duo/dx,dwo/dx) in a potential 



U d (u ,w ) 



aw. 



o 2 )/2. 



This is a conser- 



vative system, with Hamiltonian Hd = (duo/dx) /2 + 
(dwo/dx) /2 + Ud(uo,w ). For bright solitons, the field 
vanishes at infinity, and thus H d = 0. Then, as the 
functions uq(x) and u>q(x) reach their maximum values 
at x = 0, the corresponding derivatives are zero, and 
thus the necessary condition for a zero asymptotic is 
^1^=0 = Ud\ x=a = 0, which we use to relate the peak 
amplitudes: 



2 2 
U„,W m — UZ, 



0. 



(12) 



Combining Eqs. (|8|)-(|12[) 1 we obtain an approximate so- 
lution in the following simple form: 

uq{x) = u m sech. p (x / p) , wq(x) = w m sech 2 (x / p) , 



4{ Wr , 



(w m - 1) 



(w m - 1) 



(2 - w m ) 



(13) 



Here, the last relation allows us to determine w rn for an 
arbitrary a as a solution of a cubic equation, and then 
to find all other parameters as functions of a. For mis- 
matches in the interval < a < +oo, the parameter val- 
ues change monotonically in the ranges: < u m < +oo, 
1 < w m < 2, and +00 > p > 1. At a = 1 the values are 
Um = 3/v2) Wm = 3/2, p = 2, i.e. our general expression 
reduces to the exact solution (|J). Similarly, for a — > +00 
the limiting result (M) follows. 



in the linear limit: wq(x) ~ exp (— y/alrrl). Indeed, for 
a > 1, we have p < 2, i.e. the FF component is effec- 
tively wider than the SH, and thus the field decay rate 
is smaller, as is correctly predicted by Eqs. (13). In con- 
trast, for a < 1 the width of the FF component is smaller 
than that of the SH (as p < 2), and then the SH tails are 
to be determined by linear asymptotics. However, the 
solution ([l3]) overestimates the SH field localization. To 
account for this feature, we sue the soliton 'peak' from 
(U|) with a linear tail, so that the function and its first 
derivative change continuously, and obtain a more accu- 
rate expression for the SH component in the case a < 1: 

ui m sech 2 (x/p) , |x| < x a , 
w (x) = { (14) 
w (x a ) exp - x a )\ , |x| > x a , 

where x a = (p/2)\n[(2 + p^/a^ / (2 — p^/a)]. From 
Eq. (|T|) it follows that the SH profile becomes double 
scaled. That is, we predict the existence of linear tails in 
the SH component, which are effectively not trapped by 
the FF field, and carry some power: 



+ OC 



w^dx 



[w (x a )f 



a. 



(15) 



The dependence of the linear tail parameters on the nor- 
malized mismatch a is shown in Fig. 
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FIG. 1. Dependence of the SH linear tails characteristics 
on mismatch a: (a) separation from the soliton center x a and 
(b) associated normalized power P a . 



4- Soliton bound states 



3. Soliton tails 

As one of our goals is to describe the soliton tails, let us 
have a look at the far-field asymptotics that follow from 
Eqs. (|l3|). It is easy to check that the FF profile exactly 
matches the linear limit. To understand the properties of 
the SH component tails, we notice that the correspond- 
ing equation in Eqs. (|[) describes the motion of a particle 
driven by an external force Uq/2. As this expression is 
positive, the function wq(x) can't decay faster than that 



A connection between the soliton bound states and 
the linear tails has been demonstrated for many physical 
situations. For example, mutual trapping of radiating 
parametric bright solitons was shown to be possible for 
a discrete set of propagation constant values when, due 
to destructive interference, the oscillatory tails disappear 
p9| . However, in our case the localized modes are not in 
resonance with propagating linear waves |9| , and a natu- 
ral assumption is that the solitons can trap each other by 
their linear tails for a continuous range of mismatches. 
This physical picture is consistent with the results of the 
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previous numerical simulations and analytical investiga- 
tions, showing that multi-soliton states are possible only 
for a < 1 |lj,[l7],^0|,|l| . Moreover, the characteristic dis- 
tance between the neighboring solitons can be roughly 
estimated to be of order 2x a , and this expression predicts 
the non-monotonic dependence of the separation on mis- 
match a. The minimum separation should be observed 
for the mismatch corresponding to an extremum point, 
dx a /da\ — 0, and then it follows that a ma ~ 0.12, 
see Fig. [jjfa). Quite remarkably, this mismatch value 
corresponds very closely to the results of numerical cal- 
culations, see Fig. 4 in Ref. pjl. 
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of spatial solitons (a — 2). As a matter of fact, the nu- 
merical and analytical results on these plots are not dis- 
tinguishable, and that is why we show them differently, 
by continuous curves and asterisks. Plots on the right 
give corresponding relative deviations with solid lines, 
and dashed lines demonstrate the errors for the varia- 
tional solution with Gaussian profiles (SGP) [jl6 17|. We 
see that the analytical solution ([T^)-(|T^) describes the 
key soliton parameters extremely accurately. In a wide 
region of mismatch values, a > 1, the relative errors are 
smaller than 0.7% for the total power, and 0.3% for other 
characteristics. For a < 1 the deviations become larger, 
but do not exceed 3%. It is clear that our solution gives 
much more accurate results than the SGP, and the latter 
provides a slightly better estimation for the total power 
only in a narrow region of mismatches. 

At this point we would like to stress that the above- 
mentioned characteristics are not the only ones which de- 
termine the 'quality' of the approximate solutions. The 
closeness of the approximate component profiles to exact 
soliton envelopes, including the tails, is also very impor- 
tant. From Figs. ||(a) and ||(b), we see that the solution 
(p^|)-(p^) describes the profiles very accurately as well 
(left plots), unlike the SGP which matches them only 'on 
average' (see plots on the right). To characterize the ac- 
curacy numerically, we define the relative deviations as: 
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FIG. 2. Left: Comparison between numerical results (con- 
tinuous curves) and approximate analytical solutions (crosses) 
given by Eqs. (|l3|)-([L4[). The characteristics are: (a,b) maxi- 
mum amplitudes of the FF and SH, (c) total power (for the 
case a = 2), and (d) Hamiltonian. 

Right: Corresponding relative errors are shown with solid 
lines, and dotted lines present deviations for the SGP. 



5. Comparison with numerical results 



In order to determine deviations between the approxi- 
mate and exact solutions, we solved Eqs. (||) numerically. 
The corresponding dependencies of the peak amplitudes, 
total power, and Hamiltonian on the detuning parame- 
ter a are presented in Fig. || (left graphs) for the case 
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Note that these characteristics have a dB-like scale, 
i.e. smaller negative values mean better matching. In 
Fig. ||(c), the errors corresponding to solution (|l3|)-(|l4|) 
are shown on the left, and to SGP on the right. We see 
that our approach allows us to precisely describe the soli- 
ton profiles, both peaks and tails, for any mismatch a. 
That is why we were able to reveal some remarkable fea- 
tures of two-component parametric mutual trapping for 
a < 1, when the SH field configuration becomes double 
scaled [see Eq. (pd|),(|l5|) and relating discussions]. On the 
other hand, the SGP does not provide close estimations 
for the actual profiles, especially for a < 1, where the 
discrepancies increase drastically. 

Approximate solutions can be used not only in theo- 
retical studies, but also in numerical simulations. The 
unique accuracy of solution (H13)-(14) makes it an al- 
most perfect generator for 'soliton' input conditions. We 
checked that the initial propagation stage is accompanied 
by very minor oscillations, and that the associated power 
loss due to radiation is negligible [see, e.g., Fig. f|. 
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FIG. 3. Comparison of numerical results with analytical 

predictions: solution (p^)-(p^), on the left, and the SGP, on 
the right. Dotted lines correspond to the FF and solid lines to 
the SH components. (a,b) Two-wave soliton profiles at a = 4 
and a = 0.01 respectively. Continuous curves show exact nu- 
merical and crosses show approximate solutions, (c) Devia- 
tions of approximate profiles as defined in Eqs. (|lrj). 

To summarize, the approximate analytical solution, 
given in a compact explicit form by Eqs. (|l3|)-(|i"4|), de- 
scribes virtually all principal features of the soliton fam- 
ily. That is why it can be called 'almost exact'. We 
might wonder, why is it so accurate? The key point in 
the analysis was to take into account the self-similarity 
of the SH envelopes. Then, we can view Eqs. (jl3|)-(|i~4|) 
as approximate scaling transformations of the two-wave 
bright soliton family. 




FIG. 4. Almost stationary evolution of the FF (a) and 
SH (b). The initial condition is given by the approximate 
two- wave soliton solution dig) ) at a = 4, a = 2, A = l. 



B. Three- wave interaction in anisotropic medium 



Let us now investigate a more general case of three- 
wave interaction. Following Ref. [2^ ], we consider a 
double-phase matched wave interaction, such that the 
FF waves of fundamental orthogonal polarizations (i.e. 
ordinary and extraordinary) are coupled with the same 
SH component. That is, the full FF field is now vectorial, 
u = {u,v}, and the two components are determined as 



u(x, z) = U(x, z) cos [<p(x, z)] , 
v(x, z) — U (x, z) sin [(p(x, z)) , 



(17) 



where (p(x, z) is the polarization angle. For such a config- 
uration the system of normalized equations can be writ- 
ten in the following form B2j: 



.du d 2 u 

+ + u w = 0) 

oz 



.dv 



dx 2 
d 2 v 



dz dx 



j3iv + xv*w = 0, 



(18) 



dw d 2 w 



ta 



dz 



Qx2 p W + -(u 2 + v 2 ) 



Here % is the normalized component of the nonlinear 
susceptibility matrix, characterizing the v <-> w coupling 
'strength' in relation to the u <-> w parametric interac- 
tion process, and /3i is the phase mismatch between the 
orthogonally-polarized FF components. All the other pa- 
rameters have the same meaning as in Eqs. (0). 

We are interested in bright solitons supported by 
Eqs. (|l8|), and search for solutions in the form: 



u(x, z) = \u (xV\)e tXz , 
v(x, z) = Xv (xVX)e tXz , 
w(x,z) = Xwa(xV\)e 2lXz , 



(19) 



where uo and vq are the real stationary amplitudes of the 
FF waves with orthogonal polarizations, wq is the enve- 
lope of the SH component, and A > is the propagation 
constant. Then, we substitute Eqs. (|l9j) into the original 
system (18) and find a set of coupled equations for the 



soliton profiles: 

d 2 uo 
dx 2 
d 2 v 
dx 2 
d 2 w 
dx 2 



u + uqwq = 0, 



a-ivo + xvqWq = 0, 
- awo + ~(uI + Vq) = 0, 



(20) 



where the normalized mismatches are 

a = 2a + P/X, a 1 = l+f3 1 /\. 



(21) 
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Similarly to the two-wave case, the total power P and 
Hamiltonian H of system (|18|) can be found for sta- 
tionary solutions using Eqs. (|4|,|^), where Pq = P Uo + 
2aP W0 + (l/x)Pv , H = 0.4[P» + aP Wa + [a 1 / X )Pv ], 



and P Vo = /_ 



+°° „,2 



dx. 



First, let us study some general properties of Eqs. (20). 
We notice that in an isotropic medium, x — a i = 1 
and the problem reduces to the two-wave case consid- 
ered above, as the FF wave can have an arbitrary con- 
stant polarization angle p. Simple two- wave solutions are 
also possible in an anisotropic medium, but only for triv- 
ial polarizations: <p = (v Q = 0) and tp — ir/2 (u = 0). 
However, it was shown that solitons with mixed polariza- 
tions also exist p2{ . To study such three- wave paramet- 
ric coupling, we follow the same path as in the previous 
section, and, in order to understand the principal scal- 
ing properties, refer to an exact one-parameter family of 
solutions found for cti — 1/4, x = 1/3, an d a > 1 [p2[: 



u (x) = (3/V2) sech 2 (x/2) , 
vq(x) — \/3 (a — 1) sech (x/2) , 
w (x) = (3/2) sech 2 (x/2). 



(22) 



An interesting feature of this three-component solution is 
that the SH profile is the same for any a. Thus, we sup- 
pose that the SH envelope shape always remains close 
to sech 2 (x /p), just as for two-component solitons, and 
choose the trial function as in Eq. (^). The SH acts as an 
effective waveguide simultaneously for two different FF 
waves, and that is why, when solving the corresponding 
equations in Eqs. (|20|), we obtain two relations between 
the SH profile characteristics: 



1 + p" 



q(q + i) 



(23) 



2w m (a + 2p 2 ) = u 2 n + v 2 



(26) 



Then, we consider an equivalent Hamiltonian dynamic 
problem and, after matching the values Hd\ x=0 = 
-Hdlx—,+00 = 0) obtain another relation: 



,(«; m - 1) +vl l (w m - aix 1 )-ai»m = Q- (27) 



It is now convenient to turn back to the polar nota- 
tions ( p7| ) . The total FF intensity and polarization angle 
at the soliton center can be found using Eqs. (prj|),(|2T 



Ul = 2w m (a + 2p- 2 ), 

sin 2 tp m = (l - aix^ 1 ) 1 (aw^U' 



w r , 



1)- 



(28) 



Finally, approximate three-wave soliton profiles are 
given by Eqs. (|s|) , (p5| ) , with the parameters in an explicit 
form from Eqs. (|24|) , (p8|) . With no lack of generality, 
we hereafter assume that x < 1, as it is always possi- 
ble to swap the functions u and v before re-normalizing 
the physical equations from which the system ( |l8| ) origi- 
nates. Analysis reveals that three-wave solutions exist if 
the mismatches a and ot\ fulfill the following inequality: 



a (u) (ai) < a < 



+00, 



q(ai) < 1) 



(29) 



where q = y/ceip. From Eq. (|2 
find the parameters: 



it is straightforward to where 



/Oi\ 



lax 



V 



-1/2 
qa 1 



'Oil 



oti-x 



(24) 



p 2 (p - 1) ' 



M = 



4«i 



q 2 (q-iy 



(30) 



Then, the FF envelopes corresponding to these values 



u (x) = u m sech p (x/p), v (x) = v rn sech q (x / p) , (25) 



where u m and v m are the peak amplitudes. We would 
like to note that in the frames of our approach, the SH 
profile does not depend on a, as follows from Eqs. (24). 
This is quite an interesting approximate scaling property 
of the three-wave solitons, and it is exactly satisfied for 
the solution presented in Eqs. (^). 



Now we have to determine the remaining unknown pa- 
rameters, viz. the FF amplitudes u m and v m . As for the 
two-wave case, we fulfill the SH equation at x = 0: 



The corresponding polarization angles are 
tp [a -> a^] -> (i.e. v -> 0) and <p[a-> a (t,) ] -> tt/2 
(u — > 0). It is now obvious that the boundaries ( puj) 
correspond to bifurcations from a two-wave soliton to a 
three- wave one. 

An example of the parameter region (a,ai), where 
three- wave mutual trapping occurs, is shown in Fig. ^ 
for x — 1/3- We see that analytically calculated bound- 
aries given by Eqs. ( |3fj| ) (shown with dashed lines) are ex- 
tremely accurate, and almost coincide with those found 
using numerical simulations. Note that it also gives the 
correct prediction that at ol\ — ► x the two boundaries 
merge, i.e. a^ u ' v ^ — * 0. 
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FIG. 5. Shaded region shows the numerically-found ex- 
istence region for three-wave solitons; analytically calculated 
boundaries ( pS(i| ) are plotted with dashed lines. Dotted lines 
(a,b,c) show dependence of the normalized mismatches on the 
propagation constant A for f)\ = —4 and f3 — 10, —15, —25 re- 
spectively. Parameters are \ = 1/3 and a = 2. 

To understand the features of the three-wave soli- 
tons, we now recall that solutions of the original system 
( |l8| ) constitute a one-parameter family, as follows from 
Eqs. (p^). Then, according to Eq. (|2l|), a change of the 
propagation constant A corresponds to motion along a 
straight line in the parameter space (a, u-y). The limit- 
ing point for A — > +oo is (2a, 1), which always lies above 
the three-wave existence region. Thus, we find that the 
trajectory will go through this region only if (3\ < and 
(3 > 2a Pi /(l - X)- N °te that for fixed f3\ the value of 
/3 determines the inclination angle, as demonstrated in 
Fig. U by the lines (a,b). On the other hand, line (c) 
corresponds to the case when the specified conditions do 
not hold, i.e. three-wave trapping is not possible for any 
A. From this analysis it follows that the three-wave soli- 
ton always bifurcates from a two- wave state with v = 0, 
and then transforms into the other two-wave mode with 
u = 0. This process can be easily seen in Fig. [| where 
examples of the power dependence on the propagation 
constant are shown. 



3600 



2200 



(a) 
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1200 



(b) 



400 L± 
5.0 



5.3 

A 
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FIG. 6. Dependence of the total soliton power P on the 
propagation constant. Parameter values for the plots (a,b) are 
the same as in Fig. |B| for the trajectories (a,b), respectively. 
In both cases, a three-wave soliton (solid line) bifurcates from 
two- wave solutions with v = (dotted) or u — (dashed). 
Crosses show approximate analytical results. 



The plots in Fig. || demonstrate that the analytical so- 
lution gives a very precise estimate for the total power. 
It also accurately describes the soliton profiles in a wide 
region of mismatches, with an example being shown in 
Fig. |?]. A thorough qualitative comparison is, however, 
a separate task, as there exist several parameters which 
control the mutual trapping. It will be presented else- 
where. 

To summarize, we have analyzed three-wave solitons 
in an anisotropic medium. Our approach allowed us to 
predict corresponding parameter regions, and power de- 
pendence with high accuracy. Of course, some interest- 
ing aspects remain to be investigated, such as stability, 
formation of linear tails, and properties of higher-order 
modes. These problems are topics of separate study and 
thus will not be addressed here, but the results obtained 
form a background for further in-depth investigations. 



o.o I 



0.0 2.5 5.0 7.5 10.0 



FIG. 7. Three- wave soliton profiles. Left: FF components 
(uo(x) - dotted, and vo(x) - dashed). Right: SH envelope. 
Crosses give analytically-calculated amplitudes. Parameters 
are the same as in Fig. |6|(a), and A = 5.35. 



III. TWO-DIMENSIONAL SOLITONS 



Two- wave parametric spatial solitons can exist both in 
planar waveguides and in bulk nonlinear media Q . In the 
latter case, the interaction between the FF and SH waves 
is_ described by the following system of coupled equations 



du d 2 u d 2 u 

+ + ^2 + u w = 0, 
oz ox oy z 

d 2 w d 2 w 



dw 

r 

dz 



dx 2 dy 2 



(3w - 



1 



(31) 



= 0, 



where x and y are the transverse coordinates. Other 
variables and parameters are the same as for the (1+1) 
dimensional case described by Eqs. (0), with the obvious 
difference that the amplitudes depend on three coordi- 
nates, u — u(x,y,z) and w = w(x,y,z). We consider a 
spatial case close to the phase matching, i.e. a ~ 2. 
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Our goal in this section is to calculate the soliton en- 
velopes by extending the method introduced for a planar 
case. Specifically, we look for circularly-symmetric sta- 
tionary solutions of Eqs. (|3l|), which can be found with 
the help of the following substitution: 



u(x, y, z) — Xuq (rV\^ e 
(n/A) 



,2i\z 



(32) 



where r = y/ x 2 + y 2 is the radial distance in cylindrical 
coordinates, and uo(r) and wo(r) are the real normalized 
envelope functions. Then, Eqs. ( |3l"| ) reduce to 



cPuq 1 dug 

dr 2 r dr 

d 2 wo 1 dwQ 

dr 2 r dr 



u 



awQ 



u w - 
1 2 



(33) 



where the propagation constant A and the normalized 
mismatch a are introduced in the same way as in 
Eqs. (g),^). Then, the soliton power and Hamiltonian 
are expressed via the the normalized values as: 



P = \P 0l H = \ 2 (H -P ), 



(34) 



where Pq — P UQ + 2aP WQ , and 



Uq 



"+00 



f 2 f 

P Uo =2tt u a rdr, P wo = 2tt w a rdr, 

Jo Jo 

2 



Hn = 27T 



+ 00 



+ 



dug \ f dw 
dr J \ dr 
+Uq + aw — UqWq] rdr. 



(35) 



beam width. Thus, we can assume that, to some extent, 
the SH envelope form stays almost intact for a > 1, and 
the trial function is 



w Q (x) = w s F (b„r) 



(36) 



where the function F describes a characteristic SH pro- 
file, and the relative amplitude w s , together with the in- 
verse width b Sl are scaling parameters. Then, we assume 
that, similarly to the (1+1)D case ([[o|), the FF profile 
can be approximately described as 



u (x) = u s F p {b s r), 



where p is an unknown parameter. 



(37) 



To determine the SH characteristic profile, we consider 
a mismatch a = 1. This case is easier to analyze, as the 
envelopes of both mutually trapped components coincide, 
wo(r) — uo(r)/\/2 = F(r). Here the function F(r) satis- 
fies the following equation 



d 2 F 1 dF 
dr 2 r dr 



F = 0, 



(38) 



whose approximate solution was found earlier with the 
help of a Hartree-like approach |23| . However, here we 
choose to use a variational method in order to select the 
one which provides better matching. First, we present 
Eq. (pq) in a variational form: 



^=0 
SF ' 



(39) 



where S denotes the variational derivative, and L is the 
Lagrangian corresponding to the original Eq. 



Similarly to the one-dimensional case, Eqs. ( p3| ) are not 
integrable, and approximate soliton profiles were found 
using a variational method. Solutions were obtained 
for trial functions with Gaussian profiles (SGP) and de- 
scribed dependence of the Hamiltonian and power on 
mismatch a > with a reasonable accuracy The 
limitations of the SGP remain the same - rough match- 
ing of actual profiles and inadequate description of the 
soliton tails. In order to construct a better solution, we 
have to start with some approximate scaling property of 
the soliton family. The problem is that effective 'dissipa- 
tion' terms (~ in Eqs. ([33|) lead to: (i) distortion of 
soliton profiles in the center and (ii) higher far field local- 
ization, which, taken together, result in very complicated 
scaling features. This makes a general analysis quite diffi- 
cult, and thus we limit our study to mismatches a > 1. It 
has been shown that in this parameter range the FF and 
SH relative beam width changes do not exceed a factor 
of tw o fi~6l , p4| ]. On the other hand, from the structure of 
Eqs. (|33|), it follows that the balance between the second- 
order and effective dissipation linear terms, which affects 
the soliton shape, depends mainly on same characteristic 



L(F) 



+ OC 



dF 
dr 



p 2 ^ 3 



dr. 



Next, we select a trial function in the form 
F (r) = F m sech 2 (6 r), 



(40) 



(41) 



assuming that the profile is more or less close to that of 
planar solitons (||). Now we note that, according to (|3S|), 
the Lagrangian reaches a minimum at an exact solution, 
and thus, to determine the peak amplitude F m and in- 
verse width bg in the approximate expression (fi~l|), an 
extremum point of the L(Fq) integral should be found: 
dL(F )/dF m = dL{F Q )/db = 0. Solving these equa- 
tions, we obtain the parameter values: 



F„ 



15 (4 In 2- 1) 
321n2 - 11 



2.3781, 



(42) 



1 /5 41n2-l) „ oi „ 

bo = ~\ « 0.5818. 

u 2V 8n2+l 
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We made a comparison between the exact numerical and 
approximate solutions of Eq. (H^Wiven by Eqs. (ftl|),(|4^) 
and the Hartree approximation |23), and found that our 
result provides a much better matching. Thus, solution 
([i"l|) , (El) will be used in further calculation. In particu- 
lar, we derive an approximate expression for the deriva- 
tive dF/dr, which will be useful in further analysis: 



(43) 



Now, after learning some important properties of the 
characteristic SH profile, the next step is to determine 
parameters in the trial functions ( |36| , |3"7j ) . First, we sub- 
stitute these expressions into the equation for the FF 
component in system (|3|). The resulting equality can't 
be satisfied exactly for any r, but we can use Eq. 
and get approximate relations between the parameters: 



/ dF^ 








1 ^<F 2 [ 









pb 2 s [l + 46g(p-l)] =1, 

pb 2 s [l + AbKp-^F- 1 ] =w s 



(44) 



Next, we fulfill the SH equation from (|33|) at r — * 0, 
and obtain: 

w s b 2 s (l - F rn ) - aw s + ^u 2 s F%- 1 = 0. (45) 



and then it is straightforward to calculate b s , w s , and u s 
one after another, employing Eqs. ( f4^ , [45| ) . Finally, us- 
ing approximate expressions ( f4l| , [f"2| ) for the wave profiles 
at a. = 1 and scaling transformations Eqs. ( p6|j37| ), the 
solution can be written as: 



u (x) 
w (x) 



u m sech 2p (a;6), 
= WmSech 2 (xb) , 



(48) 



where u m = UsF^y/2, w m = w s F m , and b = b s b . The 
amplitudes, characteristic widths, and scaling parame- 
ter p are monotonic functions of a, and for 1 < a < 
+00 change in the limits: 1 > p > 1/2, bo < b < 
^2b 2 /{l - 2b 2 ) « 1.448, {F m V2 « 3.363) < u rn < +oo, 
F m <w m < (F m -2b$)/(l-2b 2 ) ) w 5.267. Keep in mind 
however, that, unlike the (1+1)D case, the asymptotics 
for large a are not exact. 



It is interesting to note that solution (|48[), describing 
the soliton profiles in bulk media, reduces to expressions 
( |l3| ) for the (1+1)D case if the parameters characteriz- 
ing the wave envelopes at a = 1 are chosen according to 
Eq. (§) as F m = 3/2 and b = 1/2. 

To study the accuracy of our solution, we made com- 
parisons with direct numerical calculations. The results 
are summarized in Fig. [s] in the same way as for solitons 
in planar structures (see Fig. 0). 



To determine the solution parameters, one more relation 
is needed, preferably following from a conservation law, 
as this allows us to better describe far-field asymptotics. 
The difficulty here is that the equivalent dynamical sys- 
tem described by Eqs. (33) is dissipative, and thus is not 
Hamiltonian. However, it is still possible to derive a con- 
dition similar to Eq. (^) for a planar configuration. To 
calculate an integral relation, we substitute approximate 
profiles ( p6p7| ) into the equation for the SH component 
in (|3^), multiply the equality by dF(b s r)/dr, and inte- 
grate over (0, +oo). The resulting expression contains 
an averaged dissipation term, which can be found using 
Eq. @ 



+ 00 



dr r 3 m 2 m ' 



and then the final equality is 



w s b„ 



1 



~:F„ 



1 



1 



4p + 2 



^r x =0. (46) 



Now all the parameters can be determined from the 
derived relations. First, p is found as a solution of cubic 
equation: 



8a6gp 3 + 2a (l - 6b 2 ) p 2 + [(4/3) F m - 2 
+a {Ab 2 - 1)] p + [1 - (4/3) F m ] = 0, 



(47) 
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FIG. 8. Comparison between exact numerical and approx- 
imate solutions, (^8|) and SGP, for solitons in bulk media. The 
characteristics shown are the same as in Fig. ^. 
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We also look at relative deviations, denning them in a 
way similar to ( |l6| ) as: 



+°° I, .exact ,,, a PP rox -| 2 „ 



S u = In 



S w = In 



■dr 



00 I,,, exact „, approx. i2 



rr k 



rdr 



(49) 



For larger mismatches, i.e. a > 10 1 , the deviations in- 
crease, but the errors always remain smaller than for the 
SGP. If even more accurate results are needed, it should 
be possible to start the derivation with the NLSE in the 
limit a — > +00, and then take into account terms of order 
a -1 . As for the opposite case, a < 1, linear tails may 
be expected to form, similarly to a (1+1)D configuration. 
However, such special analysis is beyond the scope of cur- 
rent paper, and we believe that these open problems will 
be addressed in further studies. 



Each corresponding dependence is shown in Fig. |9j(c) for 
our solution ( [48| ) on the left, and the SGP on the right. 
From the data presented, it follows that analytical solu- 
tion pH ) provides a very good approximation for both 
integral characteristics and soliton profiles [see also plots 
on the left in Figs, ^(a) and ^(b)]. Especially accurate 
results are observed for mismatches of the order a of 10° 
to 10 1 . This range actually corresponds to interesting 
cases from an experimental point of view, as usually soli- 
tons are observed more or less close to phase matching, 
i.e. a r~j 2a ~ 4 Jl],||]. The amazing precision is due 
to the fact that the approximate profile (^l] , ^ ) provides 
outstandingly good matching with the exact envelope, 
as shown in Fig. ^|(a). Actually, it can be claimed that 
this is an 'almost exact' solution at a = 1 for solitons 
in a bulk medium, as the deviations 5 U and S w become 
extremely small [see the left plot in Figjp(c)] . Note also 
that the approximation coefficients in [ p3f are not nearly 
as accurate as those given in ([42]). 




2 4 6 (b) 2 4 6 

r r 




FIG. 9. Comparison of numerical results with analytical 
predictions for (2+l)D soliton profiles given by: Eq. wm) on 
the left, and the SGP on the right. (a,b) Envelopes at a = 1 
and a — 10 correspondingly, (c) Deviations of approximate 
profiles as defined in Eq. (|49j). Notation is the same as in 
Fig. |. 



IV. CONCLUDING REMARKS 

We have studied the properties of two- and three- 
component quadratic bright solitons in planar waveg- 
uides and in bulk media. Very accurate, yet simple 
and compact, approximate solutions have been derived 
to describe the actual wave profiles, with a precision 
which the previously-known approximations could not 
achieve. Such amazingly good results have been obtained 
because the trial functions were chosen to correspond to 
the scaling properties of the solitons. With the help of a 
specially-developed approach, the optimal values of fit- 
ting parameters have been found in a simple explicit 
form, which wouldn't be possible with the variational 
method. 

In particular, an 'almost exact' solution has been de- 
rived which describes a whole family of two-wave planar 
solitons, accounting for all the key properties. It not only 
provides perfect estimations for integral characteristics 
(power and Hamiltonian) , but also closely matches the 
envelope profiles of both components for any mismatch. 
On the other hand, the solution allowed us to reveal the 
existence of non-oscillating linear tails, and rigorously de- 
scribe their features, which for example made it possible 
to explain and predict some peculiarities of multi-soliton 
bound states. 

We also considered three-wave coupling in an 
anisotropic medium between the orthogonally-polarized 
FF and SH components. Mismatch values, when the 
three wave solitons can exist, have been determined ana- 
lytically, and bifurcation scenarios have been described. 
Additionally, an approximate solution for soliton profiles 
has been obtained, and it provides close estimations in a 
wide parameter region. 

Quite interesting results have been obtained for the 
case of solitons in a bulk medium. Although exact an- 
alytical expressions for the two-wave profiles in a (2+1) 
dimensional case are not known, even in simpler limit- 
ing cases (e.g., in the cascading limit the resulting single 
component NLSE is not integrable), we have been able 
to derive an 'almost exact' solution for a specific mis- 
match value. General approximate expressions are also 
presented, giving very close matching in a wide range of 
detunings, covering values close to phase matching, which 
are of major interest from the experimental perspective. 
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